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A multi-scale computational model for determining plastic behavior in two-phased 
CMSX-4 Ni-based superalloys is developed on a finite element analysis (FEA) framework 
employing crystal plasticity constitutive model that can capture the microstructural scale 
stress field. The generalized method of cells (GMC) micromechanics model is used for 
homogenizing the local field quantities. At first, GMC as stand-alone is validated by 
analyzing a repeating unit cell (RUC) as a two-phased sample with 72.9% volume fraction of 
y ? -precipitate in the y-matrix phase and comparing the results with those predicted by finite 
element analysis (FEA) models incorporating the same crystal plasticity constitutive model. 
The global stress-strain behavior and the local field quantity distributions predicted by 
GMC demonstrated good agreement with FEA. High computational saving, at the expense 
of some accuracy in the components of local tensor field quantities, was obtained with GMC. 
Finally, the capability of the developed multi-scale model linking FEA and GMC to solve 
real life sized structures is demonstrated by analyzing an engine disc component and 
determining the microstructural scale details of the field quantities. 


I. Introduction 

S ingle crystal materials such as Nickel-based super alloys (e.g. ME3) are widely used in engine blades and disks 
due to their exceptional strength, fracture toughness and resistance against creep and corrosion at elevated 
temperature and other harsh environments. Nickel-based super alloys contain two-phased composite 
microstructures; y-matrix phase (predominantly an alloy of Ni, Cr and other alloying elements with a face-centered 
cubic (fee) lattice Structure) and the precipitate phase y f (a coherent, ordered Ni3Al inter-metallic of LI2 crystal 
structure). The morphology and the distribution of the y' precipitates inside y matrix influences the deformation, 
deriving excellent mechanical properties. Therefore, the volume fraction of precipitates, their mean size and spacing 
have a major role on the behavior of these alloys in different loading and temperature conditions. Comparing to the 
matrix phase, precipitates show much less nonlinear plastic behavior. As a result, the precipitates act as obstacles to 
the motion of dislocations in the matrix 1 " 3 . 

In order to predict the microstructural plastic behavior in the two -phased Ni-based super alloys, models that can 
incorporate the interaction of various microstructural features, such as grain boundaries, material phases and 
dislocations are needed. The two key characteristics associated with plasticity are the permanent strain due to 
slipping of atomic layers under shear stress and the strain hardening resulting in nonlinear stress -strain behavior. In 
the crystal plasticity models, the dislocation motion is incorporated as permanent shear strains in the various slip 
systems (slip plane and slip direction) implemented with the aid of a set of state variables 4,5 . The strain hardening is 
implemented in terms of the strength of slip systems, another set of state variables. The nonlinear evolution laws are 
defined for the permanent shear strains and strength of slip systems, along with the equilibrium equations in an 
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incremental and iterative approach. Several studies on the development of the crystal plasticity models with 
constitutive relations defining shear flow and hardening rates for Ni-based superalloys have been reported in the 
literature 6 " 14 . Studies focusing on the interaction of two phases by implementing constitutive models for the 
individual phases separately have also been reported (Refs. 6, 12 and 13). 

The major challenge in implementing microstructure based plasticity models to real sized structural components 
is the computational cost as it requires solving governing equations at both the microstructural and macrostructural 
length scales. A potential solution for addressing this challenge is the careful implementation of an appropriate and 
efficient homogenization technique to bridge between the micro-scale and the macro-scale deformation without 
significantly compromising on the accuracy of the solution, especially in the critical regions of interest. 

In this study, the generalized method of cells (GMC) 15-19 is applied as an integrated homogenization method 
which allows linking the macro-scale to the microstructure while considering the individual microstructural phases 
explicitly. To this end, a crystal plasticity model in the form of a user material subroutine is implemented in the 
constituents at the micro-scale using GMC homogenization 20, 21 . At first, the performance of the stand-alone GMC, 
in terms of solution accuracy and computational time, was evaluated by comparing with a finite element model that 
uses the same material subroutine. Finally, the multi-scale model, implemented on finite element analysis 
framework at macro-scale with element properties defined at micro-scale using GMC homogenization, was 
evaluated by analyzing a realistic engine disc component. 

II. Single-crystal Plasticity Framework 

The crystal plasticity theory is based on the assumption that any elasto -plastic deformation can be split into two 
multiplicative operations, a plastic deformation where material is deformed through the rearrangement of lattices but 
not through its own deformation, followed by an elastic deformation associated with the stretching of lattices 22-24 . 
The total deformation gradient F is then given by, 

F = F e . F p (1) 

where F v and F e are the plastic and elastic deformation gradient, respectively. 

Based on the deformation gradient definition, the total velocity gradient L is stated as 

L = F.F~ 1 =D + n (2) 

which the symmetric stretch rate D and the anti-symmetric spin tensor f2 can be decomposed into lattice and plastic 
parts (( )*, ( ) p ), respectively. 

D = D* + D p , fi = fi* + f2 p (3) 

The velocity gradient associated with the plastic deformation L p is given in terms of Schmid’s law as, 

l p = FP.pp~ 1 = Y ia y a m a ®n a ( 4 ) 

where y a is the rate of shear strain associated with the slipping of a slip system, m is the unit normal to the slip 
plane and n is the unit vector parallel to slip direction. The incremental formulation of plasticity theory is based on 
(1) the evolution of Cauchy stress on the corotational frame of reference that rotate with the crystal lattice, /* (<r), (2) 
the slipping rate y^ and (3) the strain hardening rate as given below 

n&) + o(I :D*) = C.D\ (5) 

where a is the Cauchy stress, D*is rate of stretching associated with elastic deformation and C is the tensor of elastic 
moduli. The governing equation for slipping is defined as 

Y (a) = a (a) (^) ( 6 ) 

where x a is the resolved shear stress, g a is the strength, the refers to the slipping rate when resolved shear 
stress reaches the strength. The governing equation for the strength can be modeled as 

9 {a) = £/? Kp Y {a) (7) 

where h a p is the slip hardening moduli. The sum ranges over all activated slip systems. h aa and h a p are called self 
and latent hardening moduli, respectively. 

Hardening is defined based on Pierce 25 and Asaro 26 model 

h aa = h(y ) = h 0 sec 2 h ° y (no sum on a ) ( 8 ) 

l T S -T 0 I 

Where h 0 is the initial hardening modulus, r 0 is the yield stress which equals the initial value of current strength 
5 (a) ( 0), t s is the stage-I stress (or the break-through stress where large plastic flow initiates), and y is the Taylor 
cumulative shear strain on all slip systems, i.e. 

Y = EaJoV (<r) l dt ( 9 ) 

The latent hardening modlui are given by 
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( 10 ) 


Kp = qh(r) (a * P) 

Where q is a constant. 

For the numerical implementation of single crystal plasticity model a user material subroutine in the form of UMAT 
for Abaqus 27 is used 4,5 . 


III. Generalized Method of Cells (GMC) framework 

GMC is a micromechanical formulation for predicting the overall thermo -inelastic behavior of the multi-scale 


composites. There are four steps involved in this homogenization process ’ . First the RUC should be identified 
(Figure la) and discretized. A typical RUC is consisted of N a x Np X N y rectangular subcells ( a = 1, ... , N a ; (3 = 
1, ... , Np; y = 1, ... , N y ) in the x lf x 2 and x 3 directions, respectively, as shown in Figure lb. The individual subcells 
have the dimensions of ( d a , hp, l y ) and are related to the unit cell dimensions ( d , h, l ) by, 


d = r a a =1 d a 

a) 


7 _ 7 

h — Lp =1 dp , 


l=K Y =1 dy 


( 11 ) 


b) 



N y = 2 subcells 16 

As the second step, the relationships between the macroscopic average stresses and strains with the microscopic 
ones are stated. A linearly variable displacement (i = 1 , 2, 3) is considered inside each subcell. This 

displacement profile is consisted of pure displacement components in the center of each unit cell and the micro - 
variables representing the linear dependence of the displacement to the local coordinates. 

The average strain in the RUC is defined as 

P — _}_Y N a y N P y N Y j u i ^( a Pr) 

d dhl L,a=lLl P=^ Ll Y= 1 a d 

which is based on the definition of strain tensor in each subcell 


( 12 ) 


~e[f Y) = \{.d iU ) al3Y) + diu{ apY) ) 


£ ( “ Py) as: 


-a . (13) 

Considering a general constitutive equation for the elastic-viscoplastic temperature-dependent material in each 
subcell (a, / 3 , y ), the average stress in that subcell is considered as, 

_(£*/?y) _ r (c*0y) r -(>0y) _ -/ (ajffy) _ -r(ajffy). .... 

°ij ~ L ijkl y b kl b kl b kl ) 

where the C is the elastic tensor. The eff^and are the average inelastic and thermal strain tensors, 

respectively, in each subcell. 

Based on the volumetric summation of the average stresses in all subcells, the average stress in the entire unit-cell 
can be written as, 

— 1 y Na y Nf * y Ny h h / 

°d ~ dh i^a=l^p=i^y=l damply O ( 1 ^) 

In the third step the continuity of displacement and traction is applied in an average, or integral, sense. It is 
considered that the RUCs are periodic and the interfaces of the subcells and the boundaries of the RUCs follow the 
displacement and traction continuity. Applying these conditions results in elimination of the micro -variables in the 
displacement definition and produces a set of continuity relationships. As a result, a set of N a (Np + N y + l) + 
Np(N y + l) + Ny equations stems from the strain continuity relation which can be declared as: 

A G e s =Js (16) 
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with A g and / being matrices representing the geometrical details of the subcells and those related to the unit-cell 
itself, respectively. The £ is the average strain in the unit-cell, while £ s is defined as, 

£ s = (^ 111 V..,f( JV “' JV / ?,JVy )) (17) 

The continuity of traction, would result in 

A„(e s - ej - ej) = 0 (18) 

which contains 6N a NpN y — (N a Np + N a N y + NpN y ) — (N a + Np + N y ) equations, in which A M includes the 
information regarding the material properties in each subcell. £ l s and £ T S represent the inelastic and thermal strains in 
each subcell and are defined in a similar way as Eq. 17. 

The strain and traction continuity conditions (Eqns. (16) and (17) can be combined as 


Ae s - D(e’ s + 4 ) = Ke 


where, 


(19) 


A = 






Rewriting this equation for each subcell as 

zWy) = A(“W£ + dWv\£[ + £ t s ) (20) 

allows expressing the average strain in the subcell in terms of the macrostrain, and the inelastic and thermal strains 
by defining the concentration matrices A^ a ^ and D^ Y \ 

In the fourth step, the overall macroscopic constitutive equations of the material based on the effective properties 
can be determined. The average stress in each subcell 

G = B*(£~£ ! ~£ T ) (21) 

in which the B* represents the effective elastic tensor of the composite: 

B* = zjf ! Zyh dahply ( 22 ) 

The inelastic strain tensor £ 7 and the average thermal strain tensor £ T are also obtained in a similar fashion. 

This method can be reformulated using subcell stresses rather than strains as the unknown quantities which reduces 
the computational cost significantly 28 . 


IV. Proposed mutli-scale model 

The Micromechanics Analysis Code with Generalized Method of 
Cells (MAC/GMC), developed by NASA Glenn Research Center 
that predicts the effective nonlinear response of heterogeneous 
materials using GMC based on the known behavior of the 
constituent materials is used in the present study. The 
MAC/GMC software, like most of the commercially available 
FEA packages, admits user-defined constitutive models. This 
capability was used to incorporate the single crystal plasticity 
model within the code such that the behavior of the individual 
constituent materials could be represented with this model. The 
GMC is called at finite element, integration point level to 
represent the material response of that element using NASA’s FEAMAC software ’ ’ 



Figure 2. CMSX4 two-phased Ni-based super 
alloy, (cubical precipitates can be distinguished 
from the dark matrix ) 30 


V. Numerical Implementation of two-phased single crystal materials 

In this study, CMSX-4 Nickel-based super alloy which is extensively used for turbine blades is chosen. This alloy 
has an average precipitate size of 0.45 gm and volume fraction of roughly 70% (Figure 2) (in this study 0.729). The 
composition of CMSX-4 for each phase is shown in Table l 29 . 


Table 1. Nominal composition of alloy (CMSX-4) (at %) as given in 29 



Ni 

Al 

Ti 

Ta 

Cr 

Co 

W 

Mo 

Re 

Hf 

CMSX-4 

63.1 

12.5 

1.3 

2.0 

7.6 

10.3 

2.0 

0.4 

1.0 

0.1 

Matrix (y) 

54.9 

- 

- 

- 

20 

18.6 

2.7 

0.8 

3 

- 

Inclusion(y') 

75 

23 

2 

1 

- 

- 

- 

- 

- 

- 


The inclusion is considered to be totally elastic while the matrix demonstrates elasto -plastic behavior. The elastic 
properties of individual phases 29 and the required properties for the constitutive model for the matrix 30 are also 
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shown in Table 2. The GMC results are compared with the FEA analysis to evaluate GMC as a homogenization tool 
for studying the multi-phased structures which would pave the way for enabling the multi-scale modeling for 
analyzing Ni-based superalloy engine components. 


Table 2. Nickel-Bases Super alloy (CMSX-4), elastic properties for Matrix and Inclusion 29 and Plastic properties 
for the Matrix 30 



E(G< Pa) 


G(GPa) 

n 

a 

h 0 (MPa) 

To (MPa) 

r s (MPa) 

Matrix (y) 

129 

0.391 

137 

10 

0.001 

350 

260 

775 

lnclusion(y / ) 

114 

0.405 

130 

- 

- 

- 

- 

- 


A. Numerical implementation using FEA 

Uniaxial loading with a 1 GPa load applied in the z-direction was considered. Periodic boundary conditions 
were applied in all the three directions. The geometry of the modeled two -phased single crystal (CMSX-4) is 
illustrated in Figure 2. Three different finite element meshes were considered for studying the convergence of the 
solution. 


a) b) c) 



Figure 2. Three different mesh refinements for studying the multi-phased single crystal The sample is cut in half in 
the z-direction to show the inclusion, a) one element each in matrix thickness and inclusion thickness direction 
(named as (111)) b) two elements each in matrix and inclusion thickness directions (named as (222)), and c) two 
elements in matrix thickness direction, and 9 elements in Inclusion thickness direction (named as (292)) 

Figure 2a shows a mesh containing one element in the matrix in the thickness direction while having one element for 
inclusion. This case is named as (111). The second case with two elements in each of the matrix and inclusion 
thickness direction, case (222), is illustrated in Figure 2b. In the third case, the number of elements in the inclusion 
thickness direction is increased to nine (292) (Figure 2c). 

First, the convergence of the FEA solution was studied. The average axial stress -strain curve of the samples is 
shown in Figure 3. It is observed that all the studied cases demonstrate very good agreement for both the matrix and 
inclusion, and as a result for the entire domain. In addition, in Figure 4 the axial stress distributions on the studied 
cross-section demonstrate close match. The case (111) showed uniform distribution for each component through the 
studied section, while case (292) captured the varitation; however, they are negligible relative to the average field 
distribution values. It should be noted that the stress values are maximum in the comers of the inclusion region 
which can only be observed in case (292). Therefore, for the application that requires very high accuracy in the field 
solution, the models containing higher number of elements should be utilized. 


a) 


b) 


c) 


Stress(Axial)-Strain Stress(Axial)-Strain Stress(Axial)-Strain 



strain strain strain 


Figure 3. Axial stress vs. Axial strain for the studied multi-phased single crystal for cases 111, 222, 292 in a) only 

Matrix b) only Inclusion c) whole domain 
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a) 


S,S33 
(Avg: 75%) 

_ *1 1 50c+03 

■ +1.1126+03 

■ +1 0750+03 
W- +1 038o+03 
-H +1.0000+03 

■ +9 6250+02 

■ +9 250«+02 
!■! +8.8750+02 
— U +8 500O+02 
— t +8 1250+02 
U +7 7500+02 
■- +7 3750+02 

+7 0000+02 




b) 


S, S33 
(Avg; 75%) 





Figure 4. Axial stress distribution (MPa) on the cut plane for a) (111) case b) (292) case 


B. Numerical Implementation using GMC 

Same Nickel-based super alloy (CMSX4) sample studied in the previous section is modeled using GMC based 
on the (1 1 1) mesh and (292) mesh. The results for the average axial stress-axial strain and the vonMisses stress-time 
responses, and vonMisses stress distributions on the chosen cross-section are shown in Figures 5 and 6. 


a) 
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£ 1000 
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CO 

<D 

55 500 


0. 


Stress(Axial)-Strain 


b) 


O 111-Whole 
-•-292-Whole 
o 111-Inclusion 
* 292-Inclusion 
□ 111-Matrix 
-■♦■■292-Inclusion 
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* 292-Inclusion 
. o 111-Matrix 
-♦-292-Matrix 








0.008 0.01 0.012 


"0 0.2 0.4 0.6 0.8 T 

time (s) 


Figure 5. Behavior of the studied multi-phased single crystal in GMC for cases 111 and 292 for whole domain, only 
inclusion and only matrix a) Axial stress vs. Axial strain b) vonMisses stress vs. time 
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From Figures 5 and 6, no difference is observed in the reported results from the cases (111) and (292). This was 
expected as the GMC method is a homogenization scheme, demonstrating the average value in the phase. 




X xIO - 4 X xIO - 4 

Figure 6. Results from GMC for cases 111 (Left) and 292 (Right) on same cut planes as shown in Fig. 5.2 for Axial 

stress distribution (MPa') 


C. FEA and GMC comparison 

The reported results from FEA and GMC are compared in this section to evaluate the credibility of GMC for 
studying the two-phased single crystal s like Nickel-based super allows. The results showing the axial stress against 
the axial strain and the vonMisses stress in time curves obtained from GMC and FEA (shown in Figure 7) 
demonstrate excellent agreement for both the inclusion and the matrix phases as well as the entire domain, 
a) b) 


Stress(Axial)-Strain Stress(Mises)-time 



Figure 7. Results from comparison of GMC and FEA for whole domain, only inclusion and only matrix a) Axial 

stress vs. Axial strain b) vonMisses stress vs. time 

The distributions of the individual components on the chosen cross-section are compared in Figure 8. For both GMC 
and FEA, the simplest case of (1 1 1) is chosen for the comparison. Based on the excellent agreement obtained, it can 
be concluded that the GMC is a reliable method for studying the two -phased single crystal materials. Although the 
simple geometry of this type of material allows modeling the CMSX-4 unit cell using the least required number of 
elements in FEA and GMC, still there is a significant advantage (63. 7x) in the computational cost of GMC. 
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a) 


Stress(S33) 



S, S33 
(Avg: 75%) 

S +1.150O+03 
+1.1 12e+03 
♦1 .0750+03 
+1 .0380+03 
+1.000O+03 
+9. 6250+02 
+9250o+02 
+8 8750+02 
+8 500O+02 
+8 125e+02 
♦7.750O+02 
+7.375O+02 
+7.000O+02 





S, Mises 
(Avg: 75%) 

I +1.150O+03 
+1.1 12o+03 
+1 .0750+03 
+1 0386+03 
+1 000e+03 
+9. 6250+02 
+9250O+02 
+8.8750+02 
+8.500O+02 
+8.1 250+02 
+7.750O+02 
+7.3 75o+02 
+70006+02 




c) 


r x 10 


Strain (E33) 



0.011 


0.0105 


0.01 


0.0095 


[ 0.009 




xIO 



Figure 8. Results from comparison of GMC (Left) and FEA (Right) for case 111 on the cut planes as shown in Fig. 
30 for a) Axial stress distribution (MPa') b) Axial strain (lower bound is 0.009 and upper bound is 0.011) c) 

vonMisses stress distribution (MPa) 


VI. Multi-Scale Simulation 


In order to implement a multi-scale analysis on realistic size problems, a macro-oscale FEA model is linked to 
the GMC micromechanics. The FEA determines the stress/strain on the integration points of the elements which is 
updated in each time step. These localized quantities are applied to the RUCs through GMC allowing the 
homogenization process in RUC. For the purpose of demonstrating this functionality, a segment of a real-size 
turbine disk is analyzed. The disk is 254 mm long and is exposed to a centrifugal load (p&) 2 ) with the magnitude of 

25.5 kgT< ^ d ^ (corresponding to the squared angular speed of 3e9 (^-^) 2 ). Considering the symmetry of the 

geometry and the load, 1/8 of the disk is analyzed applying appropriate symmetric boundary conditions. It is 
important to note that the symmetric assumption is valid strictly only at the macro-scale. The two-phased CMSX-4 
superalloy was used as the material model for each element in the micro-scale domain. 


8 

American Institute of Aeronautics and Astronautics 




a) 



Elem 112 


S, Mises 
(Avg: 75%) 

+ 1.229e+06 
+ 1 .1 67e+Q6 
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+9.81 2e+ 05 
+9.1 93&+D5 
+8,574e+ 05 
+ 7 . 954©+ 05 
+ 7. 335fi+ 05 
+6,71 5e+05 
+ 6 . 096 &+ 05 
+5.47G&+ 05 
+4.857e+05 


c) 
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+ 1 . 0430+06 
+9.81 2e+05 
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+7. 954 e+ 05 
+7,335e+05 
■ +6.71 5e+05 
+ 6 . 096 G+ 05 
+5.476e+05 
+4,857e+05 



tel 


I 9 MmBw<*lM8-SSl*TSi)u 





b) 




Figure 9. vonMisses stress distribution (x 10 3 MPa) in the real-size disk example a) linking macro-scale to 


micro-scale using FEA/GMC b) applying the standard FEA model with a macro-scale constitutive relationship c) 


Micro-level study at two integration points on a chosen cross-section (element number 112) 
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Figure. 9a shows the vonMisses stress distribution in the disk. To compare the macro-scale stress distribution, 
the vonMisses stress distribution obtained by using a standard macro-scale plasticity model is shown in Figure 9b. In 
the standard FEA model, the global stress-strain curve obtained from the same two-phased CMSX-4 under periodic 
boundary condition was used as the constitutive model. It is important to note that the stress/strain variation per 
element in the standard FEA follows the element interpolation while in the multi-scale model the GMC can evaluate 
the stress/strain distribution at the micro-scale consistent with the interpolated macro-scale values of the macro-scale 
finite element. 

The vonMisses stress distributions shown in Figure 9a and 9b demonstrate that the multi-scale linking, using 
GMC homogenization under finite element framework, predicts expected distribution at macro-scale. The 
vonMisses stress distributions on a chosen cross-section halfway through the depth of the unitcell for the two of the 
integration points of element 112 (shown in Figure 9a) are demonstrated in Figure 9c. The results show that the local 
field quantity distributions at the grain-level in a real-sized model can be determined using the developed multi-scale 
model which is extremely difficult with a standard FEA analysis alone due to exorbitant computational cost. 

VII. Conclusion 

A multi-scale computational model (employing FEA and GMC) that can capture the microstructural scale stress 
field for real sized problems, employing crystal plasticity constitutive relations in conjunction with the generalized 
method of cells (GMC) for homogenization of local field quantities, was developed. Two-phased materials (e.g. 
ME3) were studied by explicitly defining each individual phases. To this end, the micro-scale behavior of the two- 
phased materials were determined using GMC and compared to that obtained from standard FEA. The GMC 
analysis achieved high savings in computational cost at a minimal expense of accuracy in the components of local 
tensor field quantities. GMC showed excellent agreement with the FEA results in terms of the average stress -strain 
behavior and the field quantity distributions. Finally, the multi-scale aspect of the model, linking the microstructure 
with the two explicit phases to the macro-scale model, was demonstrated by implementing GMC as a 
homogenization tool on FEA platform, and investigating a real life turbine engine disc problem. The macro-scale 
results demonstrated similar distribution when compared to the FEA based analysis thereby verifying the method. In 
addition, it was shown that the micro-scale distribution of field quantities like vonMisses stress can be extracted on 
the cross-section of an arbitrarily selected element. The developed multi-scale model may allow engineers to tailor 
the microstructure for different structural applications. 
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